#ifndef NPDM_H
#define NPDM_H

#include <tuple>
#include <vector>
#include <set>
#include <map>
#include "npdm_typedefs.h"

//----------------------------------------------------------------------------------------------------------------------------------------------------------

class NpdmBlueprint {

  public:
    // Return sets of required operators indexed by LHS, Dot or RHS for a given DMRG sweep position
    std::map< char, std::set< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > > getSpinOpsForSite(int);
    std::map< char, std::set< std::vector< std::tuple<int,CRE_DES> > > > getSpatialOpsForSite(int);
    // Return set of irreducible spatial operator strings for a given DMRG sweep position
    std::set< std::vector< std::tuple<int,CRE_DES> > > getSpaceStringsForSite(int);
    // Split operator string into LHS, Dot and RHS blocks (overloaded on spin or spatial orbitals)
    std::map< char, std::vector< std::tuple<int,CRE_DES> > > splitIntoBlocks( int, const std::vector< std::tuple<int,CRE_DES> > & );
    std::map< char, std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > splitIntoBlocks( int, const std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > & );
    // Print functions
    void printIrreps();
    void printBlockOps();
    void printDMRGmap( );
    // Print operator string (overloaded on spin or spatial orbitals)
    void printIdx( const std::vector< std::tuple<int,CRE_DES> > & );
    void printIdx( const std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > & );

  protected:
    // Variables
    //-----------

    int kmax_;
    int pdmOrder_;
    // Set of irreducible vectors. Each vector element is one operator, stored as a 3-tuple (spatial index, spin, cre/des)
    std::set< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > idxSet_;
    // Each DMRG-ordered operator string, with original irreducible order and parity
    std::map< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > , std::pair< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> >, int > > idxMap_;
    // Bins of DMRG-ordered operator strings, according to which dot they are associated
    std::multimap< int, std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > dmrgIdxBins_;

    // Methods
    //---------

    // Combine all distinct and valid pairs of creation op strings with destruction op strings to get set of irreducible operator strings
    void buildIrrepOpStrings( std::vector< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > &, 
                              std::vector< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > & );
    // Build sets of DMRG-ordered operator strings, stored by site location
    void buildSpinDMRGmap( );
    // Get permutation parity of input operator string wrt DMRG ordered form
    int getPermParity( const std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > & ); 
    // Reorder operator string into DMRG form
    void sortIntoDMRGform( std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > & );
    // Convert spin-orbital operator string to spatial orbitals only
    std::vector< std::tuple<int,CRE_DES> > transformSpinToSpace( std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > & );

};

//----------------------------------------------------------------------------------------------------------------------------------------------------------

#endif

